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ABSTRACT 


The  development  of  the  three  dimensional,  laminar  velocity  profile 
in  the  entrance  length  of  a  rectangular  duct  is  investigated.    The  solution 
to  this  hydrodynamic  problem  is  obtained  from  the  full,  incompressible 
Navier-Stokes  equations  and  the  continuity  equation,  in  finite  difference 
form,  on  the  digital  computer  employing  the  computational  method  of 
Chorin  (1).    The  solution  yields  the  hydrodynamic  velocities  U,  V,  and 
W  and  the  friction  factor  as  a  function  of  the  distance  from  the  entrance. 
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I.     INTRODUCTION 

The  design  of  compact  heat  exchangers  for  use  in  gas  turbine 
regenerators ,  among  other  applications ,  entails  the  maximum  use  of 
the  high  heat  transfer  rates  attainable  in  the  entrance  region.    In  this 
region  the  fluid  is  undergoing  rapid  transition  from  its  uniform  profile  at 
the  entrance  to  its  fully  established  value  farther  downstream.    A  knowl- 
edge of  the  temperature  profiles  in  the  region  is  required  to  obtain  the  heat 
transfer  rates.    The  temperature  profile  develops  simultaneously  with  the 
velocity  profile  and  is  dependent  upon  it.    The  Graetz  solution  assumes  a 
fully  established  velocity  profile  and  thus  eliminates  most  of  the  complexity 
of  the  hydrodynamic  problem.    This  solution,  however,  fails  to  yield  ac- 
curate results  in  problems  involving  large  entrance  region  to  overall  length 
ratios  o 

Various  approximate  methods  have  been  devised  for  two  dimensional 
channels,  such  as  a  circular  pipe  or  parallel  plate  channels,  to  determine 
the  development  of  velocity  profiles.    Until  recently  these  solutions  have 
invariably  made  boundary  layer  assumptions.    Schlichting  (5)  obtained  a 
solution  for  the  flow  between  a  pair  of  infinite  parallel  flat  plates.    He 
used  two  asymptotic  series  solutions  ,  one  based  on  Blasiusr  solution  of 
the  boundary  layer  development  expanded  in  the  downstream  direction, 
and  the  other  based  on  the  Hagen-Poiseuille  solution  of  a  parabolic 
velocity  distribution,  expanded  in  the  upstream  direction.    He  then 
joined  these  two  solutions. 
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Recently  Wang  and  Longwell  (6)  have  employed  the  full  momentum 
equations  and  the  introduction  of  a  stream  function  to  obtain  a  solution 
for  the  velocity  development  between  parallel  flat  plates  „    The  resulting 
equations  were  cast  in  finite  difference  form  and  the  solution  effected 
numerically.    To  date,  no  exact,  three  dimensional  solutions  for  develop- 
ing velocity  profiles  have  been  obtained. 

The  purpose  of  this  analysis  is  to  present  an  exact,  finite- 
difference  solution  to  the  full  incompressible  Navier-Stokes  equations 
for  the  hydrodynamic  entrance  region  of  a  square  duct.    This  solution  is 
not  nearly  as  simple  as  approximate  methods.    Its  value,  however,  lies 
in  its  potential  accuracy  and  the  fact  that  it  is  an  exact  solution  on 
which  the  credibility  of  various  approximate  solutions  may  be  based  = 
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II.    ANALYSIS 

Ao         Governing  Equations 

The  governing  equations  are  the  momentum  (Navier-Stokes)  and  the 
continuity  equations.    The  complexity  of  the  equations  requires  the  fol- 
lowing assumptions  concerning  the  flow: 

1 .      steady  and  three  dimensional 

2        laminar 

3.      incompressible 
Through  these  assumptions,  the  energy  and  momentum  equations  become 
uncoupled ,  and  may  be  expressed  as  follows: 
Case  I: 

Momentum: 


&+^+-vfr +<»&<■ -«ff+* 


-*M    +.  d_U_+  d± 


U 


**      r&^'ym  +  t&+tf-J  (i) 

^+uib:  +  v^v;  4-um^t  = -^^+-v&+S]£+££  1 


(3) 


Continuity: 

&     <vd     ^  (4) 
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B .  Boundary  Conditions 

The  boundary  conditions  associated  with  the  equations  for  a  square 
duct  with  a  uniform  velocity  at  the  entrance  are: 
1  „      At  the  entrance  (  x  =  0) 
v  =  w  =  0 

U   =  UQ 

2„      At  the  walls  (x=0,  y  =  yw,  z  =  0,  z  =  z^; 
u  =  v  ■  w  =  0 

C .  Non-dimensionalization 

Equations  (1),  (2),  (3)  and  (4)  may  be  written  in  dimensionless  form 
through  the  introduction  of  the  following  dimensionless  variables: 


x=  a    !  "  d     ^     a  a1 

The  non-dimensional  equations  are  then: 
Momentum: 


(5) 


A2  ^  p  n  ,^ 

4T 


(6) 


<^Y     4XX     ^Y1    &* 


W  +  p  I,  ,^w  +v,^w  , .,^TL  -i£  -^i^  +^w  +  £w 

af +4  5*     S?      **J"  3?  5"?+^+3T<  (8) 


Continuity: 
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D.        Friction  Factor 

To  obtain  an  equation  for  the  friction  factor,  the  relation  between 
the  shear  stress  and  the  flow  variables  must  first  be  introduced: 

Since  the  duct  is  symmetrical,  either  derivative,  —     or  *-*-    ,  may  be 
used  in  equation  (10).    By  introducing    a  new  variable  Y|^  where  Ti  =  y  or  ^    , 
equation  (10)  may  be  written: 

The  friction  factor,  F,  may  be  related  to  the  shear  stress  by  the  equation. 

By  substituting  equation  (11)  into  (12)  the  desired  equation  for  the  friction 
factor  is  obtained: 


^=0 


By  defining  a  new  dimensionless  variable,  "VY-   y£    0  equation  (13)  may 
be  written  in  dimensionless  form 


(14) 

Y\=  O 
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III.     METHOD  OF  SOLUTION 

A.         Technique 

The  recently  formulated  technique  of  Chorin  (l)  was    adopted  for  the 
solution  of  the  governing  equations .    The  principle  of  the  method  of  solution 
suggested  by  Chorin  lies  in  the  introduction  of  an  artificial  compressibility, 
£  ,  into  the  equations  of  motion  in  such  a  way  that  the  final  results  do  not 
depend  on  £  .    The  solution  is  in  essence  a  relaxation  method.    In  employ- 
ing this  method,  it  is  necessary  to  introduce  a  modified  set  of  governing 
equations: 


Momentum: 


^X       &y        zi 


^x     ix'1    ^x   «>?■ 


:i5) 


3K[uM*w$. 


^Y     ^     $Y2    fP 


(16) 


^i      ^X1     ^Ya      ^ 


(17) 


Continuity: 


A£  4.  ii  +  AY.  t  *v£ 

^T      ^X      <^Y      ^ 


^0 


(18) 


State: 


[19) 
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In  this  set  of  auxiliary  equations  C  is  an  artificial  density, $ 
is  an  artificial  compressibility,  and  P  =  ^/&  is  an  artificial  equation 
of  state . 

Equations  (15) ,   (16),   (17)  and  (18)  may  be  put  into  finite  difference 


form  by  introducing  the  following  finite  difference  approximations: 

-U 

AX" 


&L  -  u(^ak,v,-£,t).vj(x-^><;y^j)        u(xr)  -  u(x") 
a>x  " 


AX 


|U        U(x-HAX)Y,Z,T)^-U(x-AX,y^,T)-U(x,Y,irlT^AT)-U(XlY1^JT-AT) 
^X2"  "  AXa 


.    u(xr)  +  U(X')-U(T+J-U(T') 
AX7 

The  set  of  modified  equations  put  into  finite  difference  form  and 
solved  for  the  appropriate  variables  yields  the  following  computing 
equations: 


U(T>,^k^  j~^f\ uV) - u V) 


AX1       AY        AZ' 


A.X 


^'u(Y^(y+j-.(J((Y-)v(Y-)l- 


AY 

+  DAT 
^X1 


u(?W4)-iX?M*1 


AZ  L 


U(K>U(X-)-  UCDJ  +Hfu(y^lJ(y)  .  Qd-) 


^  [u(2+)^u(r)-ur)]-f^[f(x+)-e(x^u(T-) 


(20) 
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V(T> 


1 


r 


1+  AX3      AY'5     Ai1 


-ReAT 
AX 


ru(X+)vcx^)-u(x-)v(x-) 


-3s£r 

AY  L 


vJ(vO  -v»J  -  &£  [vj(2+)v(z+)  -  w(r)v(s ) 


»  -7 


AX 


v(x^)tv(x-)-v(r))+^ 


vfr+)+v(v-)-\/fl-} 


t  |^[v^>v(r)-v(T-)]-ii  [p(Y+)-e(Y')Jtv(T") 


121) 


VJ(f)= 


1 


1+  ^*   aF^A**     I        AX    L 


ufr'MxVWMx") 


- ^fv(Y+)W(V+)-V(Y-)W(Y^ -  £^|W\*+) -Wfe") 


141 
AX' 


:at 


w(x+)+W(x-)-W(t-)|  +^r[W(Y+)+w^-)-w(r) 


^it  ,/,+< 


LT 


w(2+)+w(r)-w(T>|i  e(?+)-?fe-) 


+ 


w(r) 


(22) 


CCrt-£[uM-uoo}^ 


v(Y+)-v(Yl|-^]w(zVw(z-)fC(r) 


(2  3) 
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B.  Boundary  Conditions 

In  dimensionless  form  the  boundary  conditions  for  equations  (15)  , 
(16),  (17);  and  (18)  are: 

1,      At  the  entrance:    CX  =  0): 
V  =  W  =  0 
U  =  1 

2  At  the  walls  (Y  =  0,1     Z  =  Q',1): 
W  =  V  =  U  =  0 

3  c      At  X  =  X  entrance  length: 

W  =  V  =  0 

fully  established 

C.  Stability     . 

With  the  introduction  of  the  modified  equations  (15) ,  (16)  „  (17)  „ 
(18),  and  (19)  an  artificial  Mach  number  is  also  introduced. 

Chorin  (1)  has  stated  that ,  in  order  to  insure  stability ,  the  flow 
Mach  number  must  be  kept  less  than  one.    An  additional  requirement  that 
must  be  met  for  stability  of  the  set  is  that: 


Af^  .35G8Z  £X0)(AXorAY°r*2)S 


'/= 
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D.        Calculation  of  F 

In  order  to  calculate  the  friction  factor,  it  is  necessary  to 
introduce  a  forward  difference  approximation  as  follows: 

Using  this  approximation,  equation  (14)  may  be  written: 


(24) 


Vl=0 

Since  equation  (14)  is  applied  at  "f^  =  0 ,    U(f[)  is  identically  zero  and 
is  eliminated  from  the  finite  difference  expression, 

Equation  (24)  may  be  used  to  obtain  an  average  value  of  the  friction  " 
factor,,    An  average  value  of  -r-s?     -         may  be  determined  at  each  X-wise 
location  by  summing     U(lr\+)      at  each  grid  point  on  the  boundary  and 
dividing  by  the  number  of  points  used  in  the  summation ,    The  average 
value  of  the  slope  of  the  velocity  profile  may  now  be  used  to  calculate 

an  average  friction  factor  as  follows: 

o 
F       =    KJ  Af^  Re     £-    U"+'  <25> 

M       =    number  of  points  used 

Uwti  =    velocity  at  one  grid  point  from  the  side  wall  boundary 

The  grid  system  and  summation  points  are  shown  in  Figure  1 . 

E .         Numerical  Procedure 

A  finite  difference  solution  was  obtained,  using  the  IBM  system  360 
digital  computer.    Equations  (20),  (21),   (22),  and  (23)  involve  approximating 
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a  solution  for  the  entire  flow  field  at  time  equal  to  T  and  T     and  then 
calculating  a  new  approximation  for  time  equal  to  T    .     The  computed 
values  of  U  ,  V,  and  W  are  tested  for  continuity  in  equation  (23) ,  and 
if   C(T+)  -  £(T~)        does  not  equal  zero,  then  a  steady  state  solution 
has  not  yet  been  achieved.    T  is  then  set  equal  to  T    and  T~  equal  to  T. 
Then  the  process  is  repeated  until  a  steady  solution  is  reached.    That 
is  until       /\T  =  0       o    This  sequence  serves  to  relax  the  set  of  equations. 

Solution  also  requires  that     /^j  be  approximated  at  the  grid  points 
on  the  side  wall  boundary  since  equation  (23)  cannot  be  used  at  these 
points.    To  apply  equation  (23) ,  one  would  have  to  impose  the  condition 
that  the  flow  parameters  are  identical  on  either  side  of  the  channel  walls . 
This  assumption,  however,  is  unrealistic.    Instead,  the  following  approxi- 
mation is  introduced  at  the  side  wall  boundary: 

e(T)'=  -^v(Wf(r) 

Using  the  stability  requirements,  the  solution  to  equations  (20) ,  (21) , 
(22) ,  and  (23)  was  obtained  using  the  following  values  of  grid  parameters: 

AX=.  00 IG4        A7=.oooooi        AY=A2=.I         £=.lfc> 
The  value  of  AX   was  obtained  using  the  entrance  length  of  .0328, 
computed  by  McComas  (4) ,  for  a  hydraulic  diameter  of  1  .  0  and  a 
Reynolds  number  of  1.0.    The  fully  established  velocity  profile  at  this 
entrance  length  was  taken  from  Knudsen  (2) .    In  this  series  solution  for 
the  fully  established  value  of  U,  the  pressure  gradient     %x   was  replaced 
by  C  U    hi  I       where  F  is  the  fully  established  factor,  also  taken  from 

Knudsen  (2) . 
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The  value  of  F  computed  by  Knudsen  (56.24)  agrees  well  with  the 
value  suggested  by  Lundgren  (3)  (56 .  908)  and  McComas  (4)  (56.908). 

The  computer  program  written  in  Fortran  IV  is  included  in  the 
appendix  „ 
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III.    RESULTS 

Various  combinations  of   AX ,  AY,  A2 ,  AT,     and    £    were  used  in  the 
solution,  all  of  which  satisfied  the  stability  requirements  set  forth  by 
Chorin  (1).    It  was  also  found,  however,  that  unless  the  ratio    Lj/^y} 
was  of  the  same  order  of  magnitude  as  the  other  ratios  of  this  type ,  the 
solution  would  diverge.    The  values  of   AX,AY)A2:)  AT,      and    S     given 
previously  were  found  to  be  the  most  rapidly  convergent  values  tried. 
Profiles  were  obtained  using  these  values  and  Re  equal  to  10.    They 
are  plotted  in  Figure  2 .    The  development  of  the  centerline  velocity  is 
shown  in  Figure  3  and  Table  I,    Contrary  to  expectations,  it  was  found 
that  for  any  value  of  entrance  length,  the  flow  would  not  become  fully 

established  before  it  reached  the  X     imposed  upon  itc    This  was  true 

e 

even  though  the  X    was  much  larger  than  the  value  computed  by  McComas 
(.0328) .    In  order  to  test  the  validity  of  the  method,  a  solution  for  fully 
established  flow  in  the  duct  was  obtained  using  the  fully  established 
profile  as  boundary  conditions  at  X  =  0  and  X  =  X  entrance  length.    In 
this  case,  the  profile  at  every  X-wise  location  should  be  the  same. 
Since  the  X-wise  station    X=  Xe  A        is  the  last  station  to  be  effected 
by  the  relaxation  technique,  the  profile  at    X  =  Xe/^    is  plotted .    For  200 
iterations  the  solution  has  not  yet  converged,  and  the  probable  error  based 
on  a  calculation  of  U(Xe) /u(^%)     on  the  centerline  is  15   4  per  cent.    After 
500  iterations, the  probable  error  has  decreased  to  0.7  per  cent,  and  after 
800  iterations,  the  probable  error  is  reduced  to  0.1  per  cent.    One  can 
see  that  as  the  number  of  iterations  is  increased,  the  accuracy  of  the 
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solution  is  increased  until  a  very  high  degree  of  accuracy  (99.93  per  cent) 
is  obtained  at  1500  iterations.    This  is  shown  in  Table  II,  Figure  4  and 
Figure  5 . 

The  average  friction  factor  obtained  as  a  function  of  X  in  the  entrance 
length  is  shown  in  Figure  6  and  Table  III .    The  first  four  to  six  computed 
points  were  in  error  due  to  the  finite  value  of   AY    and    A2-      ,  which  yield 
a  finite  value  of  shear  stress  at  the  entrance .    In  actuality  the  shear  stress 
is  infinite  at  the  entrance,  due  to  the  discontinuity  in  the  velocity  at 
Y  or  i  -  0  •    In  the  limit  as  &Y\     approaches  zero,  equation  (25)  does,  in 
factff  approach  infinity.    As  the  entrance  length  is  approached,  the  value 
of  F    asymptotically  approaches  a  value  of  54.    This  value  is  in  very  good 
agreement  with  the  fully  established  value  that  was  introduced  into  the 
boundary  conditions .    It  also  is  in  good  agreement  with  the  value  suggested 
by  McComas  (4)  and  Lundgren  (3) . 
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V.    CONCLUSIONS 

As  a  result  of  this  analysis  ,  it  may  be  concluded  that: 
A  solution  using  a  finite  difference  approximation ,  and  relaxation 
type  of  scheme  suggested  by  Chorin  (l)  „  leads  to  a  convergent  and 
accurate  solution  albeit  with  substantially  more  computation  than  is 
inherent  to  an  explicit  scheme .    The  accuracy  of  the  solution  depends  on 
the  number  of  iterations  made,  the  size  of  the  finite-difference  variables 
AX  ,  &Y  -  A  2}  and  LT    ,  and  the  value  of  the  entrance  length  imposed 
upon  the  flow.    This  dependence  of  the  solution  upon  the  length  of  the 
entrance  region  is  not  what  was  expected  and  should  be  the  subject  of 
further  investigation.    The  solution  obtained  in  this  analysis  was  based 
on  a  Reynolds  number  of  1.0.    As  one  can  see,  the  computing  equations 
contain  Reynolds  number  as  a  parameter  and  therefore  depend  upon  it. 
Before  a  solution  of  this  type  may  be  used  in  the  calculation  of  heat 
transfer  in  an  entrance  region,  it  would  be  necessary  to  investigate  the 

behavior  of  this  method  at  Reynolds  numbers  that  are  more  typical  of 

3 

actual  flow  problems  (Re=  10  ) .    A  further  limitation  to  the  solution 

obtained  in  this  analysis  is  that  of  computer  time  and  storage  space. 
This  might  be  overcome  by  making  use  of  the  symmetry  of  the  duct  with 
a  suitable  modification  of  the  boundary  conditions. 
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TABLE   I 
CENTERLINE  VELOCITY 
DEVELOPMENT 


U*v, 


/U, 


0  .4815 

00328  .5314 

00656  .5816 

00984  .6321 

0   13  12  .  6  8  30 

0    16  4  0  .  7  3  43 

0   19  6  8  .  7  8  62 

02296  .83  86 

02624  8916 

02952  .9453 

0  3  2  80  1.0 
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TABLE   II 
PROBABLE  ERROR 
OF    THE 
SOLUTION 

NUMBER  OF 
ITERATIONS 

(COMPUTED)                           (ACTUAL) 

PROBABLE 
ERROR   (%) 

2    0  0 

5    0  0 

8   0  0 

1   1    0  0 

15   0  0 

.84575                                 1.0 
.99288                                 1.0 
.99908                                1.0 
.99932                                1.0 
.99933                                1.0 

15.4 
.  7 

.  0  9 
.  0    7 
.  0   7 
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TABLE    III 
CALCULATED  FRICTION  FACTOR 
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8  0 

.00164 

7   8.1 

.  0  032  8 

7  6.3 

. 0  049  2 

7  4.3 
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7  2.6 
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7  0.(7 
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6  9.05 

.01146 

6  7.3 

.0   13  12 

6  5.7 

.0  1640 

6  4.16 

.01966 

6    1.1 

.  02296 

5  8.7 

.02624 

5  5.7 

.  0278  8 

5  4.82 

.02952 

5  4.64 

.03116 

5  4.4 

.  03280 

5  4.2 
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